cache_dados <- function(chave, funcao_geradora) {
cache <- file.exists(chave)
if (!cache) {
resultado <- funcao_geradora()
saveRDS(resultado, chave)
} else {
resultado <- readRDS(chave)
}
resultado
}nH0 <- 100
nH1 <- 200
phi_parametro <- 0.2
coeficientes <- function(phi) {
# βARMA: the model from Rocha and Cribari-Neto (2009, 2017)
# is obtained by setting coefs$d = 0 and d = FALSE and error.scale = 1 (predictive scale)
list(alpha = 0, beta = NULL, phi = phi, theta = NULL, d = 0, nu = 20)
}
barma.sim <- function(n, phi, seed, y.start = NULL){
BARFIMA.sim(
n = n,
coefs = coeficientes(phi),
y.start = y.start,
error.scale = 1,
complete = F,
seed = seed
)
}
barma.phi_estimado <- function(yt, alpha = 0, nu = 20, phi = 0.1){
BARFIMA.fit(
yt = yt,
start = list(alpha = alpha, nu = nu, phi = phi),
p = 1,
d = FALSE,
error.scale = 1,
report = F
)$coefficients["phi"][[1]]
}
barma.residuos <- function(yt, phi_estimado){
BARFIMA.extract(
yt = yt,
coefs = coeficientes(phi_estimado),
llk = F,
info = F,
error.scale = 1
)$error
}# Função para forçar a execução de `BARFIMA.extract` até que ela funcione
so_quero_que_funcione <- function(func, debug = T, max = 1000, ...) {
FINALMENTE <- F
contador <- 0
while (!FINALMENTE) {
contador <- contador + 1
if (debug) print(paste("Tentativa:", contador))
resultado <- tryCatch({
func(...)
}, error = function(err) {
if (contador >= max) {
stop("Deu ruim, número máximo de tentativas atingido")
}
if (any(grepl("\\.btsr\\.extract\\(", deparse(err$trace$call)))) {
# Ignora erro de extração de resíduos
return(NULL)
}
stop(err)
})
if (!is.null(resultado)) {
FINALMENTE <- T
}
}
return(resultado)
}ewma_qcc <- function(amostra_inicial, lambda, amostra_subsequente) {
ew <- ewma(amostra_inicial, lambda = lambda, nsigmas = 3, plot = F, newdata = amostra_subsequente)
registros <- (nH0 + 1):(nH0 + nH1)
ewma <- as.numeric(ew$y[registros])
LI <- ew$limits[, 1][registros]
LS <- ew$limits[, 2][registros]
fora_de_controle <- ewma < LI | ewma > LS
total_fora_de_controle <- sum(fora_de_controle, na.rm = TRUE)
fracao_fora_de_controle <- total_fora_de_controle / length(ewma)
list(
ewma = ewma,
LI = LI,
LS = LS,
fora_de_controle = fora_de_controle,
total_fora_de_controle = total_fora_de_controle,
fracao_fora_de_controle = fracao_fora_de_controle
)
}ewma_amostras_subsequentes_gerador <- function() {
lambda <- 0.2
amostra_controle <- barma.sim(
n = nH0,
phi = phi_parametro,
seed = 1337
)
ultima_observacao <- amostra_controle[nH0]
phi_estimado <- barma.phi_estimado(amostra_controle)
residuos_controle <- barma.residuos(amostra_controle, phi_estimado)
data.frame(
novo_phi = seq(0, 0.8, by = 0.1)
) %>%
rowwise() %>%
mutate(
observacao = list(1:nH1),
dados = list(barma.sim(
n = nH1,
phi = novo_phi,
seed = 1337,
y.start = ultima_observacao
)),
residuos = list(
barma.residuos(dados, phi_estimado)
),
qcc = list(ewma_qcc(residuos_controle, lambda, residuos)),
ewma = list(qcc$ewma),
LI = list(qcc$LI),
LS = list(qcc$LS),
fora_de_controle = list(qcc$fora_de_controle),
total_fora_de_controle = qcc$total_fora_de_controle,
fracao_fora_de_controle = qcc$fracao_fora_de_controle
)
}
ewma_amostras_subsequentes <- so_quero_que_funcione(ewma_amostras_subsequentes_gerador)## [1] "Tentativa: 1"
ewma_amostras_subsequentes %>%
group_by(novo_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
.groups = "drop"
)ggplotly(
ewma_amostras_subsequentes %>%
unnest(
cols = c(novo_phi, observacao, ewma, LI, LS)
) %>%
ggplot() +
geom_line(aes(x = observacao, y = LI, color = "Limite Inferior"), linetype = "dashed") +
geom_line(aes(x = observacao, y = LS, color = "Limite Superior"), linetype = "dashed") +
geom_line(aes(x = observacao, y = ewma, color = "EWMA")) +
geom_point(data = . %>% filter(ewma < LI | ewma > LS),
aes(x = observacao, y = ewma, color = "Fora de controle"), size = 1.5, shape = 4) +
labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
scale_color_manual(
values = c(
"Limite Inferior" = adjustcolor("red", alpha.f = 0.5),
"Limite Superior" = adjustcolor("red", alpha.f = 0.5),
"EWMA" = adjustcolor("black", alpha.f = 0.8),
"Fora de controle" = adjustcolor("red", alpha.f = 0.4)
)
) +
labs(
title = "EWMA para resíduos βAR(1), com λ = 0.2",
) +
theme.base +
theme(
legend.position = "bottom",
legend.title = element_blank()
) +
facet_wrap(~novo_phi)
) %>%
plotly.baseewma_monte_carlo_gerador <- function(){
numero_de_execucoes <- 200
expand.grid(
k = 1:numero_de_execucoes,
novo_phi = seq(0.2, 0.6, by = 0.1),
lambda = seq(0.1, 0.4, by = 0.1)
) %>%
mutate(
id = row_number()
) %>%
rowwise() %>%
mutate(
amostra_controle = list(barma.sim(
n = nH0,
phi = phi_parametro,
seed = id
)),
ultima_observacao = amostra_controle[nH0],
phi_estimado = barma.phi_estimado(amostra_controle),
residuos_controle = list(
barma.residuos(amostra_controle, phi_estimado)
),
dados = list(barma.sim(
n = nH1,
phi = novo_phi,
y.start = ultima_observacao,
seed = id + 1337E4
)),
residuos = list(
barma.residuos(dados, phi_estimado)
),
qcc = list(ewma_qcc(residuos_controle, lambda, residuos)),
ewma = list(qcc$ewma),
LI = list(qcc$LI),
LS = list(qcc$LS),
fora_de_controle = list(qcc$fora_de_controle),
total_fora_de_controle = qcc$total_fora_de_controle,
fracao_fora_de_controle = qcc$fracao_fora_de_controle
)
}
ewma_monte_carlo <- cache_dados(
"cache_simulacao_ewma.RData",
\() so_quero_que_funcione(ewma_monte_carlo_gerador)
)ewma_monte_carlo_resumo <- ewma_monte_carlo %>%
group_by(lambda, novo_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
)
ewma_monte_carlo_resumoggplotly(
ewma_monte_carlo_resumo %>%
ggplot(aes(x = novo_phi, y = mean, color = factor(lambda))) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
color = "Valores de λ",
title = "Fração de pontos fora de controle por Φ e λ"
) +
theme.base
) %>%
plotly.baseggplotly(
ewma_monte_carlo %>%
ggplot(aes(x = factor(novo_phi), y = fracao_fora_de_controle, fill = factor(lambda))) +
geom_boxplot() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
fill = "Valores de λ",
title = "Fração de pontos fora de controle por Φ e λ"
) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')cumsum_qcc <- function(amostra_inicial_, se.shift, amostra_subsequente_ = NULL) {
intervalo_de_decisao <- 5
arguments <- list(
data = amostra_inicial_,
se.shift = se.shift,
decision.interval = intervalo_de_decisao,
plot = F
)
if (!is.null(amostra_subsequente_)) {
arguments$newdata <- amostra_subsequente_
}
cu <- do.call(cusum, arguments)
registros <- if (is.null(amostra_subsequente_)) {
seq(1, nH0)
} else {
seq(nH0 + 1, nH0 + nH1)
}
pos <- cu[["pos"]][registros]
neg <- cu[["neg"]][registros]
fora_de_controle_pos = pos > intervalo_de_decisao
fora_de_controle_neg = neg < -intervalo_de_decisao
total_fora_de_controle = sum(fora_de_controle_pos) + sum(fora_de_controle_neg)
fracao_fora_de_controle = total_fora_de_controle / length(registros)
list(
pos = pos,
neg = neg,
fora_de_controle_pos = fora_de_controle_pos,
fora_de_controle_neg = fora_de_controle_neg,
total_fora_de_controle = total_fora_de_controle,
fracao_fora_de_controle = fracao_fora_de_controle
)
}se.shiftQueremos que nas 100 amostras iniciais a probabilidade de um ponto fora de controle seja de 0%. Portanto,
se_shift_otimo_optimize <- function(amostra_inicial) {
for (x in seq(0, 3, by = 0.1)) {
if (cumsum_qcc(amostra_inicial, se.shift = x)$total_fora_de_controle == 0) {
return(x)
}
}
}
se_shift_otimo <- cache_dados(
"cache_se_shift_otimo.RData",
\() {
data.frame(
id = 1:100
) %>%
rowwise() %>%
mutate(
amostra_controle = list(barma.sim(
n = nH0,
phi = phi_parametro,
seed = id
)),
ultima_observacao = amostra_controle[nH0],
phi_estimado = list(barma.phi_estimado(amostra_controle)),
residuos_controle = list(barma.residuos(amostra_controle, phi_estimado)),
minimo = se_shift_otimo_optimize(residuos_controle)
)
}
)ggplotly(
se_shift_otimo %>%
ggplot(
aes(x = 0, y = minimo, group = 0, fill = "red")
) +
geom_boxplot() +
labs(
x = element_blank(),
y = "Valor ótimo de se.shift",
title = "Valor ótimo de se.shift por simulação"
) +
theme.base +
theme.no_legend
) %>%
plotly.basevalor_otimo_se_shift <- mean(se_shift_otimo$minimo)
kable(
se_shift_otimo %>%
group_by() %>%
summarise(
"Média" = mean(minimo),
"Desvio padrão" = sd(minimo)
),
caption = "Valor ótimo de se.shift",
align = 'c',
digits = 3
)| Média | Desvio padrão |
|---|---|
| 0.777 | 0.29 |
cumsum_amostras_subsequentes_gerador <- function() {
amostra_controle <- barma.sim(
n = nH0,
phi = phi_parametro,
seed = 1337
)
ultima_observacao <- amostra_controle[nH0]
phi_estimado <- barma.phi_estimado(amostra_controle)
residuos_controle <- barma.residuos(amostra_controle, phi_estimado)
data.frame(
novo_phi = seq(0, 0.8, by = 0.1)
) %>%
mutate(
id = row_number(),
) %>%
rowwise() %>%
mutate(
observacao = list(1:nH1),
dados = list(barma.sim(
n = nH1,
phi = novo_phi,
seed = 1337,
y.start = ultima_observacao
)),
residuos = list(
barma.residuos(dados, phi_estimado)
),
qcc = list(cumsum_qcc(residuos_controle, valor_otimo_se_shift, residuos)),
pos = list(qcc$pos),
neg = list(qcc$neg),
fora_de_controle_pos = list(qcc$fora_de_controle_pos),
fora_de_controle_neg = list(qcc$fora_de_controle_neg),
total_fora_de_controle = qcc$total_fora_de_controle,
fracao_fora_de_controle = qcc$fracao_fora_de_controle
)
}
cumsum_amostras_subsequentes <- so_quero_que_funcione(cumsum_amostras_subsequentes_gerador, F)
cumsum_amostras_subsequentes %>%
group_by(novo_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
.groups = "drop"
)cusum_1 <- function(){
plot <- cumsum_amostras_subsequentes %>%
filter(novo_phi <= 0.5) %>%
unnest(
cols = c(novo_phi, observacao, pos, neg, residuos, fora_de_controle_pos, fora_de_controle_neg)
) %>%
ggplot() +
geom_hline(aes(yintercept = 5, color = "Limite de decisão"), linetype = "dashed") +
geom_hline(aes(yintercept = -5, color = "Limite de decisão"), linetype = "dashed") +
geom_line(aes(x = observacao, y = pos, color = "N+")) +
geom_line(aes(x = observacao, y = neg, color = "N-")) +
geom_point(data = . %>% filter(fora_de_controle_pos),
aes(x = observacao, y = pos, color = "Fora de controle"), size = 1.5, shape = 4) +
geom_point(data = . %>% filter(fora_de_controle_neg),
aes(x = observacao, y = neg, color = "Fora de controle"), size = 1.5, shape = 4) +
labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
scale_color_manual(
values = c(
"N+" = adjustcolor("blue", alpha.f = 0.6),
"N-" = adjustcolor("darkgreen", alpha.f = 0.6),
"Fora de controle" = adjustcolor("red", alpha.f = 0.3),
"Fora de controle" = adjustcolor("red", alpha.f = 0.4),
"Limite de decisão" = adjustcolor("red", alpha.f = 0.5)
)
) +
labs(
title = "CUMSUM para resíduos βAR(1)",
) +
theme.base +
theme(
legend.position = "bottom",
legend.title = element_blank()
) +
facet_wrap(~novo_phi)
plot
}
cusum_2 <- function(){
plot <- cumsum_amostras_subsequentes %>%
filter(novo_phi > 0.5) %>%
unnest(
cols = c(novo_phi, observacao, pos, neg, residuos)
) %>%
ggplot() +
geom_hline(aes(yintercept = 5, color = "Limite de decisão"), linetype = "dashed") +
geom_hline(aes(yintercept = -5, color = "Limite de decisão"), linetype = "dashed") +
geom_line(aes(x = observacao, y = pos, color = "N+")) +
geom_line(aes(x = observacao, y = neg, color = "N-")) +
labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
scale_color_manual(
values = c(
"N+" = adjustcolor("blue", alpha.f = 0.6),
"N-" = adjustcolor("darkgreen", alpha.f = 0.6),
"Limite de decisão" = adjustcolor("red", alpha.f = 0.5)
)
) +
labs(
title = "CUMSUM para resíduos βAR(1)",
) +
theme.base +
theme(
legend.position = "bottom",
legend.title = element_blank()
) +
facet_wrap(~novo_phi)
plot
}
ggp_1 <- ggplotly(cusum_1()) %>%
plotly.base
ggp_2 <- ggplotly(cusum_2()) %>%
plotly.base
subplot(ggp_1, ggp_2, nrows = 2, margin = 0.05, heights = c(0.7, 0.3))cumsum_monte_carlo_gerador <- function(){
numero_de_execucoes <- 100
expand.grid(
k = 1:numero_de_execucoes,
novo_phi = seq(0.2, 0.6, by = 0.1),
se_shift = seq(0.1, 1.0, by = 0.1)
) %>%
mutate(
id = row_number()
) %>%
rowwise() %>%
mutate(
amostra_controle = list(barma.sim(
n = nH0,
phi = phi_parametro,
seed = id
)),
ultima_observacao = amostra_controle[nH0],
phi_estimado = barma.phi_estimado(amostra_controle),
residuos_controle = list(
barma.residuos(amostra_controle, phi_estimado)
),
dados = list(barma.sim(
n = nH1,
phi = novo_phi,
y.start = ultima_observacao,
seed = id + 1337E4
)),
residuos = list(
barma.residuos(dados, phi_estimado)
),
qcc = list(cumsum_qcc(residuos_controle, se_shift, residuos)),
pos = list(qcc$pos),
neg = list(qcc$neg),
fora_de_controle_pos = list(qcc$fora_de_controle_pos),
fora_de_controle_neg = list(qcc$fora_de_controle_neg),
total_fora_de_controle = qcc$total_fora_de_controle,
fracao_fora_de_controle = qcc$fracao_fora_de_controle
)
}
cumsum_monte_carlo <- cache_dados(
"cache_simulacao_cumsum.RData",
\() so_quero_que_funcione(cumsum_monte_carlo_gerador)
)Em média, diminuir o se.shift aumenta a sensibilidade do
CUSUM para detectar mudanças no processo.
cumsum_monte_carlo_resumo <- cumsum_monte_carlo %>%
group_by(se_shift, novo_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
)
cumsum_monte_carlo_resumoggplotly(
cumsum_monte_carlo_resumo %>%
ggplot(aes(x = novo_phi, y = mean, color = factor(se_shift))) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
color = "Valores de se.shift",
title = "Fração de pontos fora de controle por Φ e se.shift"
) +
theme.base
) %>%
plotly.baseggplotly(
cumsum_monte_carlo %>%
ggplot(aes(x = factor(novo_phi), y = fracao_fora_de_controle, fill = factor(se_shift))) +
geom_boxplot() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
fill = "Valores de se.shift",
title = "Fração de pontos fora de controle por Φ e se.shift"
) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')estatisticas_resumo_combinadas <- bind_rows(
cumsum_monte_carlo_resumo %>%
filter(se_shift > 0.7) %>%
mutate(algoritmo = "CUSUM"),
ewma_monte_carlo_resumo %>%
mutate(algoritmo = "EWMA")
)ggplotly(
estatisticas_resumo_combinadas %>%
ggplot() +
geom_line(
aes(x = novo_phi, y = mean, color = factor(lambda), linetype = algoritmo),
data = . %>% filter(algoritmo == "EWMA")
) +
geom_point(
aes(x = novo_phi, y = mean, color = factor(lambda), shape = algoritmo),
data = . %>% filter(algoritmo == "EWMA")
) +
geom_line(
aes(x = novo_phi, y = mean, color = factor(se_shift), linetype = algoritmo),
data = . %>% filter(algoritmo == "CUSUM")
) +
geom_point(
aes(x = novo_phi, y = mean, color = factor(se_shift), shape = algoritmo),
data = . %>% filter(algoritmo == "CUSUM")
) +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
linetype = "Algoritmo",
title = "Fração de pontos fora de controle",
color = "Valores dos parâmetros",
shape = "Algoritmo"
) +
theme.base
) %>%
plotly.baseestatisticas_combinadas <- bind_rows(
cumsum_monte_carlo %>%
filter(se_shift > 0.7) %>%
mutate(algoritmo = "CUSUM"),
ewma_monte_carlo %>%
mutate(algoritmo = "EWMA")
)ggplotly(
estatisticas_combinadas %>%
ggplot() +
geom_boxplot(
aes(x = factor(novo_phi), y = fracao_fora_de_controle, fill = factor(lambda), shape = algoritmo),
data = . %>% filter(algoritmo == "EWMA")
) +
geom_boxplot(
aes(x = factor(novo_phi), y = fracao_fora_de_controle, fill = factor(se_shift), shape = algoritmo),
data = . %>% filter(algoritmo == "CUSUM")
) +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
fill = "Valores dos parâmetros",
linetype = "Algoritmo",
title = "Fração de pontos fora de controle",
color = "Valores dos parâmetros"
) +
facet_wrap(~algoritmo) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')comb_monte_carlo_gerador <- function(numero_de_execucoes) {
numero_de_execucoes <- 100
expand.grid(
k = 1:numero_de_execucoes,
novo_phi = seq(0.2, 0.4, by = 0.1),
se_shift = seq(0.7, 1.1, by = 0.1),
lambda = seq(0.1, 0.4, by = 0.1)
) %>%
mutate(
id = row_number(),
"(se;λ)" = factor(paste0(se_shift, ";", lambda))
) %>%
rowwise() %>%
mutate(
amostra_controle = list(barma.sim(n = nH0, phi = phi_parametro, seed = id)),
ultima_observacao = amostra_controle[nH0],
phi_estimado = barma.phi_estimado(amostra_controle),
residuos_controle = list(barma.residuos(amostra_controle, phi_estimado)),
dados = list(barma.sim(n = nH1, phi = novo_phi, y.start = ultima_observacao, seed = id + 1337E4)),
residuos = list(barma.residuos(dados, phi_estimado)),
ewma = list(ewma_qcc(residuos_controle, lambda, residuos)),
cumsum = list(cumsum_qcc(residuos_controle, se_shift, residuos)),
fora_de_controle = list(
ewma$fora_de_controle & (cumsum$fora_de_controle_pos | cumsum$fora_de_controle_neg)
),
fracao_fora_de_controle = sum(fora_de_controle) / nH1
)
}
comb_monte_carlo <- cache_dados(
"cache_simulacao_comb.RData",
\() so_quero_que_funcione(comb_monte_carlo_gerador)
)comb_monte_carlo_resumo <- comb_monte_carlo %>%
group_by(`(se;λ)`, novo_phi, se_shift, lambda) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
)
comb_monte_carlo_resumoggplotly(
comb_monte_carlo_resumo %>%
ggplot(aes(x = novo_phi, y = mean, color = `(se;λ)`)) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
color = "Valores dos parâmetros (se;λ)",
title = "Fração de pontos fora de controle por Φ"
) +
theme.base
) %>%
plotly.basemelhores <- c("1;0.3")
ggplotly(
comb_monte_carlo_resumo %>%
filter(`(se;λ)` %in% melhores) %>%
ggplot(aes(x = novo_phi, y = mean, color = `(se;λ)`)) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φλ",
y = "Fração de pontos fora de controle",
color = "Valores dos parâmetros (se;λ)",
title = "Fração de pontos fora de controle por Φ"
) +
theme.base
) %>%
plotly.baseggplotly(
comb_monte_carlo %>%
filter(`(se;λ)` %in% melhores) %>%
ggplot(aes(x = factor(novo_phi), y = fracao_fora_de_controle, fill = `(se;λ)`)) +
geom_boxplot() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
fill = "Valores de `(se;λ)`",
title = "Fração de pontos fora de controle por Φ e se.shift"
) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')